//
// Created by 11311 on 2022/12/25.
//

#ifndef HEURISTIC_ACO_HPP
#define HEURISTIC_ACO_HPP
#include<iostream>
#include <cmath>
#include <numeric>
#include <algorithm>

using namespace  std;
namespace ag {
    typedef struct DICT{
        int key;
        double val;
    }dict;

    int roll(dict *first, dict *last){
        // sums
        double val = accumulate(first, last, 0.0, [](double a, dict &b)->double {return a + b.val;});
        // sort
        sort(first, last, [](dict &a, dict &b)->bool{return a.val<b.val;});
        dict *p = first;
        // random
        double pr = (double)rand() / RAND_MAX;
        double sums = 0.0;
        //core
        while (1){
            if(p == last)break;
            sums += p->val / val;
            if(sums > pr || abs(sums - pr) < 1e-6)return p->key;
            p++;
        }
        printf("pr = %.4f and val = %.4f and sums = %.4f\n", pr, val, sums);
        cout<<"roll function has been filed!\n";
        return -1;
    }

    void memwarn(){
        cout<<"Out of Memory!\n";
        exit(-1);
    }
    template<typename T>
    T *create_arr(int size){
        T*res = (T *)malloc(size * sizeof(T));
        if(res == nullptr)memwarn();
        for(int i = 0; i < size; i++)res[i] = (T)0;
        return res;
    }

    template<typename T>
    T** create_arr(int r, int c){
        T** res = (T**)malloc(r * sizeof(T*));
        if(res== nullptr)memwarn();

        int T_block_size = sizeof(T);
        for(int i = 0;i < r; i++){
            T *p = (T *)malloc(c * T_block_size);
            if(p == nullptr)memwarn();
            for(int j = 0; j < c; j++)p[j] = (T)0;
            res[i] = p;
        }
        return  res;
    }

    template<typename T>
    T** RandData(int r, int c, int low, int high){
        T** res = create_arr<T>(r, c);
        for(int i = 0; i < r; i++){
            for(int j = 0;j < c; j++){
                res[i][j] = (T)(rand() % (high - low) + low);
            }
        }
        return res;
    }

    class aco {
    private:
        int ant_num;
        double alpha;
        double beta;
        double rou;
        int node_num;
        int now;
private:
        void print_log(int, int,double, double);
        void one_scale() const;
        double best = 1e+16;
    public:
        double **infos;
        double **graph_path;
        double **one_scale_path;

        int **individuals;
        double *fitness;

    public:
        aco(int, int, double, double, double, double**);
        ~aco();

    public:
        double cal_fitness(int *);
        void init();
        int unit_select(int, int);
        void select();
        void update();

        double run();

    public:
        double &operator()(int epoch);
    };

    void aco::one_scale() const {
        double val = -1;
        for(int i = 0; i < node_num; i++){
            for(int j = 0; j < node_num; j++){
                if(graph_path[i][j] > val)val = graph_path[i][j];
            }
        }

        for(int i = 0; i < node_num; i++){
            for(int j = 0; j < node_num; j++){
                one_scale_path[i][j] = graph_path[i][j] / val;
            }
        }

    }
    aco::aco(int ant, int node, double a, double b, double r, double **graph) {
        ant_num = ant;
        node_num = node;
        alpha = a;
        beta = b;
        rou = r;
        graph_path = graph;
        now = 0;

        one_scale_path = create_arr<double>(node_num, node_num);

        individuals = create_arr<int>(ant_num, node_num);
        fitness = create_arr<double>(ant_num);
        infos = create_arr<double>(node, node);

    }
    aco::~aco() {
        for(int i = 0; i < ant_num; i++)free(individuals[i]);
        free(individuals);
        free(fitness);
        cout<<"Memory Has been Cleaned\n";
    }
    double aco::cal_fitness(int *c) {
        double val = 0.0;

        for(int i = 0; i < node_num; i++)val += graph_path[c[i]][c[(i+1) % node_num]];

        return val;

    }

    int aco::unit_select(int idx, int end) {
        int buff[node_num];
        for(int i = 0;i < node_num; i++)buff[i] = 0;
        for(int i = 0; i < end; i++)buff[individuals[idx][i]] = 1;
        int count = 0;

        dict values[node_num - end];
        for(int i = 0;i < node_num; i++){
            if(buff[i] == 1)continue;
            int last = individuals[idx][end - 1];
            values[count].val = pow(infos[last][i], alpha) * pow(1.0/one_scale_path[last][i], beta);
            values[count].key = i;
            count++;
        }

        int node = roll(values, values+count);

        return node;
    }

    void aco::select() {
        for(int i = 0; i < ant_num; i++){
            int node = unit_select(i, now);
            individuals[i][now] = node;
        }

    }

    void aco::update() {
        double **buffs = create_arr<double>(node_num, node_num);

        for(int i = 0; i < ant_num; i++){
            double val = cal_fitness(individuals[i]);
            fitness[i] = val;
            for(int j = 0; j < node_num; j++){
                int now = individuals[i][j];
                int next = individuals[i][(j+1) % node_num];
                buffs[now][next] += one_scale_path[now][next] / val;
            }
        }

        for(int i = 0; i < node_num; i++){
            for(int j = 0; j < node_num; j++){
                infos[i][j] = (1 - rou) * infos[i][j] + buffs[i][j];
            }
        }

    }
    void aco::init() {
        double sums = 0.0;
        for(int i = 0; i < node_num; i++)sums += one_scale_path[i][(i+1)%node_num];
        for(int i = 0; i < node_num; i++){
            for(int j = 0; j < node_num; j++){
                if(i == j)infos[i][j] = 0.0;
                else{
                    infos[i][j] = (double)ant_num / sums;
                }
            }
        }


        for(int i = 0; i < ant_num; i++){
            individuals[i][0] = i % node_num;
        }
        now++;

    }

    double aco::run() {
        for(int i = 0; i < node_num - 1; i++){
            select();
            now++;
        }
        update();
        now = 1;

        double val = 1e+16;
        for(int i = 0;i < ant_num; i++){
            if(fitness[i] < val){
                val = fitness[i];
            }
        }

        return val;
    }

    double & aco::operator()(int epoch) {
        // initialize some relative parameters
        one_scale();
        init();

        // clock
        clock_t st, et;
        st = clock();

        // running
        for(int i = 0; i < epoch; i++){
            double val = run();
            et = clock();
            if(val < best){
                double t = (double)(et - st) / CLOCKS_PER_SEC;
                print_log(i+1, epoch, val, t);
                best = val;
            }
        }

        et = clock();
        print_log(epoch, epoch, best, (double)(et - st) / CLOCKS_PER_SEC);

        return best;
    }

    void aco::print_log(int step, int epoch, double val, double t) {
        printf("[%5d/%5d] ---- UB = %6.2f\t TotalTime = %6.4f\n", step, epoch, val, t);
    }
} // ag

#endif //HEURISTIC_ACO_HPP
